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ABSTRACT 

Context. Explicit numerical computations of hypersonic or super-fast differentially rotating disks are subject to the time-step constraint 
imposed by the Courant condition, according to which waves cannot travel more than a fraction of a cell during a single time-step 
update. When the bulk orbital velocity largely exceeds any other wave speed (e.g., sound or Alfven), as computed in the rest frame, 
the time step is considerably reduced and an unusually large number of steps may be necessary to complete the computation. 
Aims. We present a robust numerical scheme to overcome the Courant limitation by improving and extending the algorithm previously 
known as FARGO (Fast Advection in Rotating Gaseous Objects) to the equations of magnetohydrodynamics (MHD) using a more 
general formalism. The proposed scheme conserves total angular momentum and energy to machine precision and works in Cartesian, 
cylindrical, or spherical coordinates. The algorithm has been implemented in the next release of the PLUTO code for astrophysical 
gasdynamics and is suitable for local or global simulations of accretion or proto-planetary disk models. 

Methods. By decomposing the total velocity into an average azimuthal contribution and a residual term, the algorithm approaches 
the solution of the MHD equations through two separate steps corresponding to a linear transport operator in the direction of orbital 
motion and a standard nonlinear solver applied to the MHD equations written in terms of the residual velocity. Since the former step is 
not subject to any stability restriction, the Courant condition is computed only in terms of the residual velocity, leading to substantially 
larger time steps. The magnetic field is advanced in time using the constrained transport method in order to fulfill the divergence-free 
condition. Furthermore, conservation of total energy and angular momentum is enforced at the discrete level by properly expressing 
the source terms in terms of upwind Godunov fluxes available during the standard solver. 

Results. Our results show that applications of the proposed orbital-advection scheme to problems of astrophysical relevance provides, 
at reduced numerical cost, equally accurate and less dissipative results than standard time-marching schemes. 

Key words. Methods: numerical - Accretion, accretion disks - Protoplanetary disks - magnetohydrodynamics (MHD) - Turbulence 



1. Introduction 

The physics of accretion flows has received a great deal of at- 
tention over the past few decades, and several numerical inves- 
tigations have contributed to broaden our current understanding 
of the problem in its diverse aspects. Representative astrophysi- 
cal scenarios involve magnetically driven turbulence in accretion 
disks around compact objects, which are typically modeled via 
a local approach (the so-called shearing -box model of |Hawley] 



et al. 1995 i or global disk simulations, that have recently be- 



come more appealing owing to increases in the available com- 
putational power see, for instance, (see, for instance Beckwith 
|et al.|201 l[[Flock et al.|2"0TT]|Sorathia et al.|2012| >. Likewise, in 
the context of planet formation, simulations of proto-planetary 
disks have become an important tool for the study of the dynam- 
ics of the gas and its impact on either particles or planets (see 
|Kley et al.|2009l|Uribe et al.|201 1| and references therein). 
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A common ingredient of all models is a highly supersonic 
or superfast, rotating, sheared flow. In standard explicit numeri- 
cal calculations, the time step is restricted, for stability reasons, 
by the Courant-Friedrichs-Lewy (CFL) condition ( Courant et al. 
|1928[ ), which is formally stated as 

Al 

Af~C fl — , (1) 

^max 

where Al is the cell length, A max is the fastest signal speed, and 
C„ - the Courant number - is a limiting factor typically of < 1. In 
essence, the CFL condition given in Eq. ([T| prevents any wave 
from traveling more than a fraction of a grid cell. This leads to 
a drastic reduction in the time step, whenever the orbital speed 
becomes considerably greater than either the sound or fast mag- 
netosonic velocities, and results in excessively long computa- 
tions. A typical occurrence is encountered in accretion or proto- 
planetary disks where the Keplerian velocity is the dominating 
flow speed. 

In the context of hydrodynamics, |Masset| ( |2000| l presented 
the first fast Eulerian-transport algorithm for differentially rotat- 
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ing disks (FARGO) in which the timestep is limited by the veloc- 
ity fluctuations around the mean orbital motion rather than by the 
total azimuthal velocity itself. The algorithm is based on the sim- 
ple recognition that if the orbital velocity is written as the sum 
of a constant average term plus a residual, the evolution operator 
can be implemented as the composition of a linear transport op- 
erator and a nonlinear step involving the residual velocity. Since 
the former step yields a uniform shift along the orbital direction, 
the CFL condition is determined during the nonlinear step by the 
magnitude of the residual velocity. The algorithm may be viewed 
as a temporary change to the local co-rotating frames placed at 
different radial positions of the differentially rotating disk. The 
advantage of using such a strategy is twofold: on the one hand, 
it reduces numerical dissipation as the Keplerian flow is now 
treated as a mean flow in the equations and, on the other hand, it 
allows larger time steps to be taken thus reducing the computa- 
tional costs (as long as the residual velocity remains subsonic). 

The ideas and concepts behind FARGO were extended to 
magnetohydrodynamics (MHD) by Johnson et al. (2008b) and 
Stone & Gardiner (2010) in the context of shearing-box model 
which is a local approximation describing a small disk patch em- 
bodied by a Cartesian box co-rotating with the disk at some fidu- 
cial radius. Because of the difficulties inherent in simulating an 
entire disk, local shearing-box models have largely contributed 
to much of what is presently known about the magneto-rotational 
instability and its implication for the transport of angular mo- 
mentum in disks. Nevertheless, several authors have stressed, 
over the past decade, the importance of global disk models, 
specially in the context of magneto-rotational instability, e.g., 
Armitage| (fT998>, |Hawley| ([2000^ , |Fromang & Nelson] ( [200%^ 
and Regev & Umurhan ( 2008). Owing to the increased numeri- 
cal power, high-resolution simulations of global disk models are 
now becoming amenable to investigation. In this context, a first 
implementation of the FARGO scheme for MHD was presented 
by|Sorathia et al. ( 2012| l as part of the Athena code ( |Stone et aL 
2008) > in cylindrical geometry. 

Here we propose, using a somewhat different and more gen- 
eral formalism, an extension of the original FARGO scheme 
to the equations of MHD for different systems of coordinates 
in the context of Godunov-type schemes. The algorithm differs 
from the original one in that it uses a fully unsplit integration 
scheme when solving the MHD equations in the residual veloc- 
ity. Source terms arising from the velocity decomposition are 
treated carefully in order to restore the conservation of total en- 
ergy and total angular momentum. The magnetic field is evolved 
via the constrained-transport (CT) scheme, which preserves the 
divergence-free condition to machine accuracy at all times. The 
proposed FARGO-MHD scheme is implemented in the PLUTO 
code for astrophysical gasdynamics (Mign one et al.|2007 2012| l 
and will be available with the next code release (4.0). The algo- 
rithm has been fully parallelized so that domain decomposition 
can be performed in all three coordinate directions, thus allow- 
ing an efficient use of a large number of processors. 

The paper is organized as follows. In section |2j we briefly 
review the basic MHD equations and the time-step limitations 
imposed by standard explicit time-stepping methods (j p.l| i. The 
new FARGO-MHD scheme is presented in 5 2.2 and consists 
of i) a standard nonlinear step (§2.3 1 formulated so as to pre- 



serves total energy and total angular-momentum conservation 
(£2.3 1 and ii) a linear transport step (52.4i where fluid quan- 
tities are simply advected in the direction of orbital motion. 
Numerical benchmarks and astrophysical applications are pre- 
sented in Section[3] while conclusions are drawn in Section|4] 



2. Description of the FARGO-MHD scheme 

2.1. Standard approach and limitations 

We begin by considering the ideal MHD equations written as a 
nonlinear system of conservation laws 



dp 
dt 



+ V ■ (pv) = 



d(pv) 

dt 



+ V ■ [pvv T - BB T ] + Vp, = -pVd> 



dE 

— + V • UE + p,)v -(v-B)B] = -pv ■ V«D 

dt 



dB 

dt 



Vx(vxB) = 0, 



(2) 



(3) 



(4) 



(5) 



where p is the mass density, v is the fluid velocity, O is the grav- 
itational potential, B is the magnetic field, p, = p + B 2 /2 is the 
total pressure accounting for thermal (p), and magnetic (B 2 /2) 
contributions. The total energy density E is given by 



1 2 1 „2 

-l + 2 PV+ 2 B 



(6) 



where an ideal equation of state (EoS) with specific heat ratio 
F has been used. Dissipative effects are neglected for the sake 
of simplicity, although they can be easily incorporated in this 
framework. 

In the usual finite-volume approach, Eq. Q to Eq. <|5j are 
discretized on a computational mesh spanned by the grid indices 
i, j, and k corresponding to the location of a given cell in a three- 
dimensional (3D) coordinate system. We label the generic unit 
vector along a given axis with e c { and consider, in what follows, 
Cartesian ({£</} = x,y,z), cylindrical polar ({<;,/} = R,<f>,z), and 
spherical ({ed} = f, 0, <f>) coordinates. Individual cells have spa- 
tial extents given by Ald,ijk, where d labels the direction. 

The time step Af is determined by the CFL condition in Eq. 
([TJ, the precise form of which can depend on the employed 
time-stepping scheme. In PLUTO, one can take advantage of 



either the corner-transport upwind (CTU) method of Colella 
( | 1990) 1 and |Gardiner & Stone| ( |2008| > or resort to a Runge-Kutta 



(RK) discretization, where the spatial gradients are treated as the 
right-hand side of an ordinary differential equation (method of 
lines). Both algorithms are dimensionally unsplit but stable un- 
der somewhat different Courant conditions. Here we consider the 
6-solve CTU and the second-order Runge-Kutta (RK2) scheme, 
which have the same number of Riemann problems per cell per 
step. If A^dim = 2, 3 is the number of spatial dimensions, one has 
( |Beckers| 19921 |Toro| 1999D 



C fl = Af max 



max 

d 



Af 

C a = — — max 

^dim y* 



z 

V d 



|Vrfl + Cj4 

Aid 

\Vd\ +Cf4 

AL 



1 



A^dim - 1 

1 



(CTU) 



(RK2) 



(7) 



dim 



where Vd and c/.d are the fluid velocity and fast magnetosonic 
speed in the direction d. For highly super-fast, grid-aligned 
flows, Eq. (|7]l shows that both CTU and RK2 yield comparable 
time steps in two-dimensions, while, in three-dimensions, RK2 
is expected to give time steps that are approximately twice as 
large. 

Nevertheless, numerical computations of advection- 
dominated flows for which Vd » c/ may result in small time 
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steps yielding unusually long computations and eventually 
suffering from an excess of numerical dissipation. This kind 
of scenario worsens in the case of a Keplerian flow in polar or 
spherical coordinates, since the cell length becomes increasingly 
smaller towards the inner radius where the orbital flow velocity 
is faster. 



2.2. An orbital advection scheme for MHD 

To overcome the limitations outlined in the previous section, 
we now assume that the fluid velocity may be decomposed as 
V = v' + W, where w is a solenoidal velocity field and v' is the 
fluctuation or residual. Eqns. (|2]i through |5]) may then be written 

as 



dp 



d(pv') 
dt 

dE' 



+ V • (pv') + w ■ Vp = (8) 
+ V • (pvV - BB + \p t ) + w ■ V {pv') = S' m - pVO (9) 



at 



+ v- 



(£' + p,)v' - (v' • B)B 



+ W-VE' = 



S' E (10) 



8B 

dt 



V x (v' x B) - V x (w x B) 



-pv' ■ VcD 

0, (11) 



where I is the identity matrix and E' is the residual energy den- 
sity, defined as 



E' = 



P 1 , 1 1 1 

T-l 2 H 2 



(12) 



and the two additional source terms in the momentum and en- 
ergy equations may be written, after some algebra, as 



and 

S' E = -pv' ■ (v • Vw) + B ■ (B ■ Vw) . 



(13) 



(14) 



Explicit expressions for these terms in different systems of coor- 
dinates can be found in Appendix § jA| 



We note that, Eqns (|8]l through (Hi have a general validity 
since no assumption has been made about the magnitude of the 
residual velocity V . However, neither the residual momentum 
nor the residual energy density are conserved quantities because 
of the appearance of the additional source terms S' m and S ' E given 
by Eq. ( 13 i and ( 14 1. Since nothing has of course changed at the 



continuous level, both the total energy and momentum must still 
be conserved. The situation is different, however, at the discrete 
level where differential vector identities are satisfied only within 
the truncation error of the scheme. This is discussed in more 
detail in § T2~3~T1 

If we look at the individual components of the system of 
Equations dH) through (Hi, each evolution equation has the form 



dq 
~dt 



+ V- 



F q +w 



(15) 



where q e (p,pv' ,E' ,B), F q is the usual MHD flux written in 
terms of the residual velocity v'. The source term S q appearing 
on the right-hand side of Eq. (15i accounts for several contri- 
butions that include gravity, the additional term S' m (or S E ), and, 
in the case of curvilinear coordinate systems, geometrical factors 



arising from differential operators. The last term on the left-hand 
side of Eq. ( 15 1 describes the linear transport of q with advection 



speed w. If w coincides with the average azimuthal velocity, this 
term has the simple effect of pushing fluid elements along their 
orbit. 

An effective approach to solve Eq. (15i is to use operator 



splitting and separate the solution into a standard nonlinear step 
that does not include the linear advection term w ■ Vq 



(16) 



and a linear transport step corresponding to the solution of 



Uq) ■ 



dq 

~dl 



+ W ■ Vq = . 



(17) 



The two operators are separately described in §2.3 and §2.4 re- 
spectively. 

When the orbital motion coincides with one of the coordinate 
axes, e.g. w = w ■ y, the exact solution of ( 17 1 can be formally 



expressed as q(y, t" + ) = q(y - wAt",t"), which corresponds to 
a non-integer translation of a row of computational zones. Since 
this operation is decomposed into an integer shift plus a remain- 
der, the linear transport step is unconstrained and can be carried 
out for arbitrarily large time steps. 

The time step is thus ultimately determined by the Courant 
condition during the nonlinear step, that is, by solving the regular 
MHD equations (plus additional source terms) written in terms 
of the residual velocity v'. Provided that Iv'J + Cf j <s; \w\, the 
CFL condition in Eq. |7]) now results in appreciably larger time 
steps, since the dominant background orbital contribution has 
now been removed from the computation. 



2.3. Standard nonlinear step 

During the standard nonlinear step, we solve Eqns ([8]) through 



(Hi without the linear advection operator (w ■ V). The resulting 
system is equivalent to the regular MHD equations written in 
terms of the residual velocity v' - v - w plus the additional 
source terms S' m and S' E given, respectively, by Eq. ( 13 i and Eq. 
•0- 

This entitles us to employ the formalism already developed 
for the solution of the MHD equations using any stable conser- 
vative finite-volume or finite-difference discretization available 
with the PLUTO code. Accordingly, we update each conserved 
quantity q using the general building block 



q" 



At" 



V-F q + S q 



(18) 



which is the discrete version of Eq. ( 16 1. In Eq. ( 18 1, the flux 



F q follows from the solution of Riemann problems at zone in- 
terfaces, while the divergence term is expressed as the sum of 
two-point difference operators in each direction d. Dropping the 
index q for simplicity and considering a generic flux F with com- 
ponents Fd — F ■ (id, we use 



Z 



1 

A 7 ^ 



(A d , + F d , + - A d ,-F d ,-) 



(19) 



where F d j± is the flux through the upper (+) and lower (-) cell 
boundaries with surface normal Cd, while Ad and 'Vd are the in- 
terface areas and cell volume. 

The magnetic field is evolved using the constrained transport 
formalism (see, e.g., the papers by Mignone et al. 2007; Flock 
|etal.|2010p . 
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2.3.1. Conservative formulation 



As anticipated in 5 2.2 the source term S q contains both point- 



wise contributions (e.g. gravitational or curvilinear terms) and 
terms involving the derivatives of the orbital velocity w. Point- 
wise source terms not involving spatial derivatives are added to 
the right-hand side of the equations in a standard explicit way. 
Conversely, source terms contributing to both the azimuthal mo- 
mentum component and the energy equation require some addi- 
tional considerations. To this end, we consider the net change in 
the total azimuthal momentum during a single time-step update. 
Using, for simplicity, Cartesian coordinates and neglecting grav- 
ity, a straightforward combination of the equations of continuity 
and the y-component of momentum yields 



(p Vy y +l = ( pVy ) n - At 



V • F m >, + wV ■ F p + pv • Vw 



(20) 



Although the sum of the last two terms on the right-hand side is 
equal to V ■ (wF p ) at the continuous level, this may not necessar- 
ily hold in a discrete sense. At the numerical level one should in- 
deed expect \wV -F p+pv' -Vw-V -{wF p )\ = 0(Ax s ), where 0(Ax s ) 
is the truncation error of the scheme. As a consequence, total 
(linear or angular) momentum and, similarly, total energy den- 
sity will be conserved only at the truncation level of the scheme. 

We, instead, seek a numerical discretization that allows to re- 
store the exact conservation of total linear (for Cartesian geom- 
etry) or angular (for polar grids) momentum and energy, also at 
the numerical level. For this purpose, we step back in the deriva- 
tion from Eqns ( |2"|)-( |5|l to ([8]l-( 1 1 1 and note that the source term, 
as given in Eq. \l3) , of the momentum equation simply results 
from the algebraic manipulation of 



which shows that, by adding the product of Eq. ( f23] l and w to 
Eq. (24 1, one obtains the conservative update of the total mo- 



mentum with corresponding flux F„ h + wF p . Since similar ar- 
guments hold by combining the density and momentum equa- 
tions with Eq. ( 25 1, the resulting discretization enforces the con- 



servations of both total momentum and total energy in both a 
continuous and a numerical sense. The corresponding expres- 
sions for polar and spherical coordinates are a straightforward 
extension of these equations leading to both total angular mo- 
mentum and total energy conservation. We report them in jj |A.2| 
and §A.3| respectively. We point out that the importance of a 



conservative formulation of the equations has already been rec- 
ognized by other authors, e.g. |Kley| ( |1998) >, who has shown that 
the inclusion of non-intertial forces as a source term may lead to 
erroneous results. 



2.4. Linear transport step 

As anticipated, the solution of the linear transport equation in 
Eq. 



17 1 consists of a uniform shift |w Af| of the conserved vari- 



able profiles in the direction of orbital motion. If the flow is pre- 
dominantly aligned with one of the coordinate axes, say y, we 
can write the solution of Eq. ( 17 1 for t e [t",t" +l ] as q(y, t) - 
q(y - w(t - t"), t"). Extensions to polar and spherical coordinates 
are straightforward, provided that y — > (p and wt — > Of. 
Integration of Eq. 



17 1 for a time step At and over a cell size 



Ay gives 



1 

Ay 



ry J+ i ry h i 

1 q n (g)dt- 2 q"(&d{ 

Jy . i -wAt \)y . i -wAt 



(26) 



S' m = -V ■ (pv'w) + w>V ■ (pv') - pw ■ Vw , 



(21) 



where the last term is identically zero in Cartesian coordinates 
and only contributes to the radial-momentum source term in 
polar and spherical coordinates. Similarly, one can show, after 



some algebra, that the source term in Eq. ( 14 1 of the energy equa- 
tion results from 



S' E = *> 



V • (pv'v' -BB+ pv'w) 



w 



V • (pv') 



pw 



v' + w ■ (pv'v' - BB) 



(22) 



where q"(g) = q(i;, f) and w does not depend on y. The inte- 



grals on the right-hand side of Eq. ( 26 1 can be converted into 
a finite sum corresponding to an integer shift of cells m - 
floor(wAf/Ay+ 1/2) plus a fractional remainder 6y = wAt-mAy. 
The final result may be cast as a two-point flux difference scheme 



where j„ 



Ay 
j - m and 



i p 

J+5 Sy 1. 



0' 



(27) 



(28) 



This suggests that Eq. pT[ ) and (22i may be conveniently ex- 
pressed in terms of the density and y-momentum upwind fluxes 
already at disposal during the conservative update. For instance, 
we update density, (residual) azimuthal-momentum, and (resid- 
ual) energy as 



p «+l _ p n 

At' 1 



V ■ F r 



(m;)" +1 - ( m ;r 

At n 

(£")" +1 - (E'T 
At" 



■ V ■ (F„ lv + wF p ) + wV ■ F f 



-V-lF £ + wF fflr + yF p | + 
+wV • (F fflv + wF p ) - ^-V • F f 



(23) 



(24) 



(25) 



is the upwind numerical flux. We remark that the finite summa- 
tions corresponding to an integer cell shift (implicitly contained 
in the integral of Eq. [26| do not need to be explicitly computed 
since most terms cancel out when taking their difference. We 
also note that, since \6y\ < Ay by construction, the scheme is 
always stable regardless of the choice of At. 



The fluxes given by Eq. ( 28 1 may be computed to an arbitrary 
order of accuracy by assuming a piecewise polynomial represen- 
tation of the data inside the cell. If, for instance, we assume a 
piecewise linear distribution of q inside the zone, then the inte- 
gral in Eq. ( 28 1 takes the form 



Aqj I 6y 



A gj+i /, ^ Sy 
q j +1 -—[ 1 + Ay- 



if 6y > , 
if 6y < , 



(29) 



4 



Mignone et al: A FARGO Scheme for MHD 



where Aqj is computed using a standard slope limiter. Equation 
(29 1 corresponds to the well-known second-order MUSCL- 
Hancock scheme of van Leer| ( |1984) >. Similarly, using a piece- 
wise parabolic distribution inside the cell we have 



with H, + i computed from Eq. <J28J> with q 



-B z . In the previ- 



H 



9M 



9j+i,- 



Ay 



$9j ^ x2 



6y\ 
Ay) 



Ay 



■ <5 2 9j+ 



Sy_ 
Ay 



if Sy>0, 



if 6y<0, 



ous equations, we have dropped, for simplicity, the full subscript 
notations when redundant and kept only the half-increment no- 
tation to represent the different magnetic and electric field com- 
ponents. 

We point out that the updating formulas for B x and B z could 
^gbe implemented in exactly the same way as done for Eq 



where q^ ± are third (or higher) order rightmost and leftmost in- 
terpolated values with respect to the cell center, 6qj = q^+ - qj- 
and 6 2 qj = q^+ - 2qj + q^. This is essentially the PPM scheme 
for a linear advection equation ( Colell a~& Woodward| 1984 1 and 
is our default choice, unless stated otherwise. 

2.4.1. Magnetic field transport 

In the constrained transport formalism, the three components of 
magnetic field are discretized on a staggered mesh and evolved 
as surface averages placed at the different zone faces to which 
they are orthogonal. In this sense, we locate the components of 



B by means of a half-integer subscript, that is, B x , 



jjk> B y,i,j+l,k 



and B„ , j k+ i to denote the magnetic field components centered 
on the x, y, and z faces, respectively. 

The evolution is carried out using a discrete version of 
Stokes' theorem applied to Eq. (Ill, where a line integral of the 



electric field is properly evaluated at zone edges. This guaran- 
tees that the divergence-free condition of the magnetic field is 
also maintained to machine precision during the linear transport 
step. During the transport step (i.e. when v' - 0), this leads to 
the following discretization of Eq. (JTTJ 



8 Z 



B 



,n+] 



B 



= B" . , 



Z,k+k 



-& z 



Ayj 



a x 



J+2 



'+5J+3 



Axi 



(31) 



(32) 



px 



j- 



Ay,- 



(33) 



where B x = — f wB z dt and & z = J wB x dt are the x and z com- 
ponents of the time-integrated electromotive force computed, in 
analogy with the results of the previous section, as 



'+5J+5 



+ < 



Ay £ B" 




B" 



if m > 



if 



if m < 



(34) 



with H: + i computed from Eq. d28]i using q — B x and 



j+ 



,k+k 



+ < 



B" 



f=j„M 



B" 



if m > 



if m = (35) 



for m < 



27 1 



since the differences of fluxes along the y-direction, defined by 
Eq. ( 3 1 1 and ( 33 1, leads to the cancellation of most terms leaving 
only the upstream values. On the other hand, cancellation does 
not occur when updating B x , since Eq. ( 32 1 involves differences 



of terms along the x and z direction corresponding to the winding 
of the field under the action of the shear. For this reason, the 
full summation must now be retained in Eq. ( [34| and ( [35] ) when 
evaluating the time-integrated electromotive force. 

For the sake of completeness, we also report the correspond- 
ing update expressions for other geometries. In polar coordinates 
(R, <p, z), we have 



r,n+l 



B n+1 



B 'rM 



= B" 



-Ri 



J+kM 5 



Az t 



it j+ 



- A', & 



(36) 



(37) 



Jt+J 



ARi 



pR 



A<f>j 



(38) 



where & R = -jQ.B z dt, 8 Z = jQB R dt (with Q. = w/R) are 
computed similarly to Eq. ( [34} and ( [35} . 

Likewise, in spherical coordinates (r, 8, <p) we have 



B" 



6,j+\ 



B" . , 

r,i+ j 



D/l 



■Ml 



i+k,k-i 



A4> k 



Acf> k 



sin# 



- sin 6: i£ 



AO; 



(39) 



(40) 



(41) 



+ sin( 



where EI = fO.B g dt,& e = - 
computed similarly to Eqs. (34 



riArt 



Q.B, dt, and O = w/(r sin 0) are 



and ( 35 1 



2.5. Expected speedup 

In general, the expected time-step gain is problem-dependent but 
can also be affected by several factors such as the magnitude of 
the orbital velocity, the grid geometry, and the cell aspect ra- 
tio. An estimate of the speedup may be directly computed from 
Eq. |7]) by taking the ratio of the time step obtained with orbital 
advection to that obtained without. Specializing, for example, 
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to polar coordinates and RK2 time-stepping, one can estimate a 
time-step increase 



max 

At F ijk 



At, 



\vr\ 



AR 



c f , R \v'^ + w\+c u \v : \+c Lz 

+ +— + 



RA(f> 



Az 



max 

ijk 



\VR\+ Cf,R 



AR 



K\ + C U Iv z | + c u 
+ 



RA(p 



Az 



(42) 



If the characteristic signal fluctuations are of comparable mag- 



nitude, |v r | + Cf , 



c f,<P 



\v z \ +c fy 



A' and are such that 



A! <c \w\, the previous expression simplifies to 



max 

AtF >Jk 



1 M + l 1 

AR RAcf> Az 



At s 



max 

ijk 



1 11 
AR R~A~4> Az 



(43) 



where M = |w|/|/l'| should be of the same order as the fast 
magnetosonic Mach number. The previous estimate shows that 
a sensible choice of the grid resolution has a direct impact on 
the expected gain: a finer resolution in the azimuthal direction 
(RA(f> <K AR) favors a larger gain, whereas cells with a smaller 
radial extent may become less advantageous. 

2.6. Parallel implementation 

In parallel computations, the simulation domain is decomposed 
into smaller sub-grids, which are then solved simultaneously 
by several processing units. For time-explicit calculations, only 
boundary data stored in the ghost zones have to be exchanged 
between neighboring processors. With FARGO-MHD, however, 
this approach presents some difficulties since the linear transport 
step may involve shifts of an arbitrary number of zones along the 
orbital direction. This implies that data values could in princi- 
ple be exchanged between all processors lying on the same row 
of a parallel domain decomposition, thus involving many more 
communications than a regular exchange of ghost zones between 
adjacent processors. 

In the tests and applications presented here, we observed that 
the maximum zone shift m m . dx hardly ever exceeds the typical 
grid size N$ of a single processor. As efficient parallel appli- 
cations require N$ > 8 - 16, we can safely assume that the 
condition m max < N# is virtually always respected and does 
not represent a stringent requirement for most applications. In 
such a way, parallel communications are performed as a cyclic 
shift between neighboring processes only along the orbital di- 
rection and the computational overhead becomes approximately 
1 + m max /« g times larger than a regular boundary call (where n g 
is the number of ghost zones) in any direction. Notice also that, 
when the orbital velocity does not change sign across the do- 
main, the amount of subtask communication can be halved since 
information always travels in the same direction and data values 
need to be transferred from one processor to the next following 
the same pattern. 

3. Numerical benchmarks 



We present a number of hydrodynamical and magnetohydrody- 
namical test problems where the proposed orbital advection al- 
gorithm is directly compared with the standard traditional inte- 
gration method. Unless stated differently, we make use of the 
ideal EoS with F = 5/3 and the PPM method, given by Eq. < |30) , 
is the default interpolation during the linear transport step. 
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Fig. 1. Local potential vorticity for the hydrodynamical vortex test after 
20 orbits at the resolution of 1024 x 4096. Top panel: standard integra- 
tion. Bottom panel: orbital advection. 

3.1. Vortex dynamics in a two-dimensional Keplerian flow 

We begin by considering the dynamical evolution of an anti- 
cyclonic vortex in a two-dimensional (2D) Keplerian disk using 
polar coordinates {R, (f>). The initial background state is defined 
by constant density and pressure equal respectively to p = 1 
and p = l/(rM 2 ), where M - 10 is the Mach number at 
R = 1, The disk rotates with angular velocity Q.(R) - 7? -3 ^ 2 
under the influence of a point-mass gravity <t> = -l/R and 
fills the computational domain defined by 0.4 < R < 2 and 
< <p < 2n. With these definitions, the disk scale-height is given 
by H(R) = c s /Q.(R), where c s = 1 jM. A circular vortex is super- 
posed on the mean Keplerian flow and initially defined by 



6v R 



/cexp 



x 2 +y 2 



' cos (f) sin <p 
- sin d> cos d> 



-y 

x l 



(44) 



where h = H(l)/2 is the size of the vortex in terms of the lo- 
cal scale-height at R = 1, k = — lis the vortex amplitude, and 
x — R cos (f> - Ro cos (po and y = R sin (f> - Rq sin <po are the 
Cartesian coordinates measured from the center of the vortex 
initially located at Rq = 1 and <po = n/4. 

We solve the hydrodynamical equations (no magnetic field) 
using the second-order Runga-Kutta time-stepping with C a - 
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Fig. 2. Time evolution of the vortensity minimum (top, normalized to 
its initial value) and relative variation in the total integrated vorticity 
(a> z (t)) = {z ■ V x v(/)> (bottom) for the hydrodynamical vortex prob- 
lem. Dotted, dashed, and solid lines refer to computations obtained at 
the resolutions of 256 x 1024, 512 x 2048, and 1024 x 4096 grid zones, 
respectively. Black and red colors correspond to the the standard inte- 
gration algorithm and FARGO-MHD, respectively. 



0.4 using linear reconstruction and follow the system until t - 
2Q0n, i.e. for 100 revolutions of the ring at R — 1. Computations 
are carried out with and without orbital advection at three differ- 
ent resolutions corresponding to 256 x 1024 (low), 512 x 2048 
(medium), and 1024 x 4096 (high), yielding approximately 
square cells in the proximity of the vortex. To enforce conserva- 
tion, we ensure that no net flux is established across the domain 
and therefore choose reflective boundary conditions at the inner 
and outer radial boundaries and impose periodicity along the <p 
direction. 



Table 1. Average time step for the hydrodynamical vortex test problem 
at different resolutions using the standard scheme (2 nd column) and with 
orbital advection (3 rd column). The time step gain is reported in the last 
column. 



NrXN$ 


A? (Standard) 


At (FARGO-MHD) 


Gain 


256 x 1024 


1.15 10~ 3 


1.40 10^ 2 


12.07 


512 x 2048 


5.73 10~ 4 


6.61 10- 3 


11.53 


1024 x 4096 


2.84 10" 4 


3.13 10- 3 


10.99 



After a few revolutions, the vortex experiences a nonlinear 
adjustment to its final persistent structure. This process is ac- 
companied by the emission of spiral density waves that radiate 
away the energy excess of the initial state ( Bodo et al. [2007 



as £ = [V x (v - QR<l>)\ z lp after 20 revolutions using the standard 
integration algorithm (top) and the proposed FARGO scheme 
(bottom) at the highest resolution. The two snapshots are in good 
agreement, although the simulation using orbital advection re- 
quired only ~ 37,000 steps compared to the ~ 444,000 steps 
required by the standard computation. We found indeed that the 
time step increased, on average, by a factor between ~ 1 1 and 
~ 12 with orbital advection, as reported in Table [T] 

As an indicative measure of the dissipation properties of the 
scheme, we compare, in the top panel of Fig [2] the time his- 
tory of the local vortensity minimum for the selected grid sizes. 
An increase in the grid resolution favors a slower vortex decay 
leading to the formation of stable, long-lived vortex structures 
( Bodo et al.|2007 1. The same effect is obtained, at an equivalent 
resolution, by employing orbital advection. 

Similarly, we inspected the relative variation in the total in- 
tegrated vorticity presented as a function of time in the bottom 
panel of Fig [2] Strictly speaking, we note that vorticity in a 2D 
compressible flow is a conserved quantity only if the fluid is 
barotropic, i.e., if p — p(p). Nevertheless, our results indicate 
that the amount of generated vorticity decreases from * 18% (at 
the lowest resolution) to « 10% (at the largest) using the stan- 
dard integration scheme (black lines), while it remains smaller 
than 1% when orbital advection is employed (red lines). 

Lastly, we verified that both total energy (including the grav- 
itational contribution) and angular momentum are conserved to 
machine accuracy, as expected from the conservative formula- 
tion of our scheme. 



3.2. Magnetohydrodynamical vortex in a shear flow 

In the next example, we investigate the stretching of a 2D MHD 
vortex in a super-fast sheared flow using Cartesian coordinates. 
The initial condition is similar to the one used by Mig none et al.| 
(2010) and consists of a uniform-density (p = 1), background 

Mtanh(x)/2, where M is 



shear flow with velocity profile v 
the sonic Mach number. The magnetic field and pressure distri- 
butions are, respectively, given by 



(B x ,B y ) = (-y,x)fie' 



(l-r 2 )/2 



1 1 

P= f + 2' 



t-^,,2 



p £ (l-r l ), (45) 



Figure[T]shows the local potential vorticity, or vortensity, defined 



where r 2 = x 2 + y 2 and, for the present example, we adopt p = 
0.02 and T = 5/3. 

We choose the square region x,y e [-5,5] as our computa- 
tional domain and impose periodic boundary conditions in the y- 
direction while applying reflecting conditions in the x-directions. 
These choices ensure that no net flux is established across the 
domain boundaries so that density, momentum and energy are 
globally conserved. We note that, in the absence of shear (for 
M - 0), the previous configuration is an exact solution of the 
ideal MHD equations, representing a circular magnetic field- 
loop supported against a pressure gradient. In the present con- 
text, however, we choose M = 20 and follow the stretching of 
the initial loop configuration as time advances. Computations are 
carried on 512 2 grid zones using the PPM scheme with Courant 
number C a = 0.8. 

Results with and without the FARGO-MHD algorithm are 
shown in Fig. [3] where we compare the magnetic pressure distri- 
bution for t = 0, 1 , 2, 3. The two solutions look indistinguishable, 
as also visible from a one-dimensional cut of the y-component of 
magnetic field at t — 3, shown in the left panel of Fig [4] A plot of 
the time step is shown in the right panel of the same figure, where 
one can see that orbital advection attains a larger gain (> 9) dur- 
ing the first four revolutions and decreases to * 4.5 towards the 
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Fig. 3. Evolution of the magnetic pressure distribution, P,„ = (B\ + Sj)/2 (xlOO), for the MHD vortex test computed using the standard MHD 
update (top row) and with FARGO-MHD (bottom row) using 5 12 2 zones. From left to right, we show the evolution from time t = to t = 3. 



end. In terms of CPU time, the standard integration took approx- 
imately 2 hours and 31 minutes whereas only 26 minutes were 
required using the orbital advection scheme. This gives, on av- 
erage, a speedup of * 5.8. 

To check whether total energy and angular momentum are 
conserved, we repeated the integration for a much longer time, 
corresponding to 100 revolutions, and added a random horizon- 
tal velocity component of the form 



,-W2) 2 



(46) 



where % is a random number between and 1. We compare, 
in Fig. [5] the conservative and the non-conservative variants of 
the orbital advection algorithm by plotting the volume-average 
of the y component of momentum and the relative variation in 
the total energy as time advances. As expected, the former con- 
serves momentum and energy to machine accuracy, while the 
latter exhibits increasing deviations from zero as the computa- 
tion proceeds. 

3.3. Shearing-box models 

The shearing-box approximation ( Hawley et al.|1995 i provides 
a local Cartesian model of a differentially rotating disk. By ne- 
glecting curvature terms, one focuses on the evolution of a small 
rectangular portion of the disk in a frame of reference co-rotating 
with the disk at some fiducial radius where the orbital frequency 
is Qq. In this approximation, the orbital motion is described by a 
linear velocity shear of the form w = -qQ.Qxy, where 



1 



1 ,ilogQ 2 (r) 

2 d log r 



(47) 



is a local measure of the differential rotation. 

The large-scale motion of the disk is described by assuming 
that identical boxes slide relative to the computational domain in 
the radial direction, a requisite implemented by the shearing-box 
boundary condition, which enforce sheared periodicity. For any 
flow quantity q not containing v v , this is formally expressed by 



q(x, y, z) -> q(x + L x , y - qD^L x t, z) 

v y (x, y, z) -> v y (x + L x ,y- qQ®L x t, z) + q£loL x , 



(48) 
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Fig. 4. Top panel: horizontal cuts at y = of the y component of mag- 
netic field at / = 3 for the MHD vortex test in Cartesian coordinates 
using the standard integration method (solid line) and FARGO-MHD 
(plus signs). Bottom panel: time-step variation during the computation. 
Solid and dashed lines correspond to standard and FARGO-MHD com- 
putations, respectively. 
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Fig. 5. Temporal evolution of the volume-averaged y component of 
momentum (black) and fractional average energy variation (green) for 
the MHD vortex with random perturbation using the non-conservative 
(dashed lines) and conservative (solid line) versions of the orbital ad- 
vection scheme. 

where L x is the radial (x) extent of the domain. The implemen- 
tation of the shearing-box boundary conditions is done similarly 
to Gressel & Ziegler (2007) and employ the same techniques 



described in this paper. 

We report, for the sake of completeness, the basic equations 
in Appendix A.l.l| and refer th e reader to|Hawley et al.|(|1995| l, 
Gressel & Zieglerf p007] l and |Stone & Gardiner| ( [20 10| i for a 
thorough discussion. In the tests below, we employ an isothermal 
EoS, assume a Keplerian flow (q = 3/2), and adopt the unsplit 
CTU+PPM scheme with a Courant number C a = 0.45. 

3.3.1. Compressible MHD shearing waves 

As a first benchmark, we consider the evolution of compress- 
ible, magnetohydrodynamical shearing-waves in an isothermal 
medium using the shearing-box model. These waves can be re- 
garded as the extension of fast and slow modes to a differen- 
tially rotating medium and provide the natural basis for wave 
decomposition in the linear theory of rotating MHD shear flows 
Pohnson|2007l l. 



We follow the same configuration used by Johnson et al. 
( 2008b]) and |Stone & Ga rdiner (2010) and consider a Cartesian 
domain of size [-1/4, 1/4] 3 with N x N/2 X N/2 grid zones. The 
initial conditions, the details of which may be found more pre- 



cisely in Johnson et al. (2008a), consists of a background equi- 
librium state about which a plane wave perturbation is imposed 



p = 1 + dp cos (k ■ x) 



w + 



-dp cos (k ■ x) 



(49) 



A = (o,0,^-^J + 6A(2,-l,5)sin(k-x) 



where w = -q£loxy is the background shear flow, H - c s /Q.q is 
the scale height, and k - An/H{-2, 1, 1) is the initial wavenum- 
ber, whereas 



dp - e 




(50) 




1 .0 1 .5 2.0 
t (orbits) 

Fig. 6. Azimuthal magnetic-field perturbation as a function of time for 
the compressible MHD shearing-wave problem. The exact solution is 
plotted as a solid black line and the dashed lines correspond to compu- 
tations obtained at the resolutions of A'; = 8 (blue), 16 (green), and 32 
(red) zones in the vertical direction. 



are the perturbation amplitudes of density and vector potential 
(e = 10~ 6 ). Magnetic field is initialized from B — V x A. We set 
the isothermal sound speed to c s = \, while the orbital frequency 
is fio = 1 ■ 

In Fig [6] we compare, at different resolutions N - %, 16, 32, 
the temporal evolution of the azimuthal field perturbation 



6B,(t) = 2 j [B y (t) - B y (f)] cos (k(t) ■ x) dV , 



(51) 



where k(t) — k + qQ.k y tx is the time-dependent wavenumber and 
By(f) = c s (l/5 - gQof/lO) is the analytical expectation, kindly 
provided to us by G. Mamatsashvili. Our results are in excel- 
lent agreement with those of | Johnson et al.| (|2008b]> and [Stone 



& Gardiner ( 20 1 0) ) showing t 



lat at the highest resolution of 32 
zones per wavelength, the solution has converged to the semi- 
analytical prediction. 

3.3.2. Nonlinear MRI 

In the next example, we investigate the nonlinear evolution of the 
magneto-rotational instability (MRI) in the shearing-box model. 
The initial condition consists of a uniform orbital motion in the 
y direction w = -qQ.Qxy of constant density p — 1 and zero 
net-flux magnetic field in the vertical direction 



B 



sin(2^)i 



(52) 



where B = 400 and c s is the isothermal sound speed. Following 
[ Gressel & Zieglerl ( |2007] i, |Johnson etaL] ( |2008b| l, and |Stone & 
|Gardiner| ( |2010| , we set c s = Qq = 10 so that the disk scale- 
height is H = Qo/cj = 1 and one orbital period corresponds to 
2jt/D.q in code units. The size of the computational domain, in 
units of H is [-4,4] x [-4,4] x [-5, 5] and 32 zones per scale- 
height are employed. We assume periodic boundary conditions 
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Fig. 7. Density distribution in the shearing-box computation after 40 
rotation periods. Top and bottom panel refers to the model with and 
without FARGO-MHD, respectively. 



in the azimuthal (y) and vertical (2) directions, while shifted pe- 
riodicity is imposed at the radial boundary. We carry out two 
sets of computations, with and without orbital advection, at the 
resolution of 256 x 256 x 32 zones. 

The evolution is characterized by an initial transient phase 
where all the relevant physical quantities grow rapidly within a 
few orbits. Subsequently, the flow dynamics becomes nonlinear 
and the system settles down into a saturated regime accompa- 
nied by the formation of typical trailing spiral density waves as 
shown in Fig.|7]at t - (i.e. after 40 revolutions). We note 

that, owing to the turbulent chaotic behavior of the MRI, a direct 
comparison between these structures is not really helpful but that 
one should instead inspect spatially averaged physical quantities. 
Fig. [8] shows the temporal evolution of the volume-integrated 
magnetic energy density (top panel) and Maxwell stresses (mid- 
dle) obtained with and without FARGO-MHD for the first 40 
orbits. The two computations are in excellent agreement with 
the absolute value of the (normalized) Maxwell stresses reach- 
ing, in both cases, the approximate constant value w 0.01. The 
time step is plotted in the bottom panel of Fig [8] where, after 
the initial transient phase, it is shown that the employment of or- 
bital advection results in a value * 4.3 larger than the standard 
calculation, in accordance with the results of |Stone & Gardiner] 
< |2010) 
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Fig. 8. Evolution of the MRI in the shearing-box model during the first 
40 orbits. From top to bottom, the three panels show, respectively, the 
volume-integrated magnetic energy-density, the Maxwell stresses, and 
the time increment used in the computations without (solid line) and 
with FARGO-MHD (dashed line). 



3.4. Disk-planet interaction 

We simulate the interaction of a planet embedded in a viscous 



global disk as in Kley et al. (20091. The 3D (r,8,<p) computa 



tional domain consists of a complete annulus of the protoplan- 
etary disk centered on the star, extending from r m ; n = 0.4 to 
'"max = 2.5 in units of ro = flj up = 5.2 AU. In the vertical direc- 
tion, the annulus extends from the disk's midplane (at 6 - 90°) 
to about 7° (or 6 = 83°) above the midplane. The mass of the 
central star is one solar mass M, = M Q and the total disk mass 
inside [r mm ,r max ] is Mdisk = O.O1M . For the present study, we 
use a constant kinematic viscosity coefficient with a value of 
v = 10 15 cm 2 /s, which corresponds to an c-value of a — 0.004 
at ro for a disk aspect ratio of H — 0.05 r sin 9. The resolution 
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Fig. 9. Top: time step as a function of orbits for the disk- 
planet interaction problem using the standard scheme (black) 
and FARGO (red). Bottom: change in angular momentum as a 
function of time. 



of our simulations is (N r , No, A^) = (256, 32, 768). At the radial 
boundaries, the angular velocity is set to the Keplerian values, 
while we apply reflective, radial boundary-conditions for the re- 
maining variables. In the azimuthal direction we use periodic 
boundary conditions hold while zero-gradient is imposed at the 
vertical boundaries. 

The models are calculated with a locally isothermal configu- 
ration where the temperature is constant on cylinders and has the 
profile T(R) oc R l with the cylindrical radius R = rsind. This 
yields a constant ratio of the disk's vertical height H to the radius 
R, which was set to 0.05. The initial vertical density stratification 
is approximately given by a Gaussian 



p(r, 0) = p (r) exp 



(n/2 - Of r 2 



2H 2 



(53) 



Here, the density in the midplane is po(r) oc f~ , which leads 
to a S(r) oc r~ l l 2 profile of the vertically integrated surface den- 
sity. The vertical and radial velocities v g and v r are initialized to 
zero, while the initial azimuthal velocity v v is given by the equi- 
librium between gravity, centrifugal acceleration, and the radial 
pressure gradient. The simulations are performed in a frame of 
reference co-rotating with the planet using a second-order Runge 
Kutta scheme with linear reconstruction and Courant number 



C fl = 0.25. Following Kley (T998JI, we adopt a conservative 
treatment of the Coriolis force in the angular momentum equa- 
tion. 



1 



Fig. 10. Logarithmic density map for the disk-planet interaction 
after 100 orbits in the equatorial plane using FARGO (top) and 
the standard integration scheme (bottom). 

The planet is described by adding a term to the total gravita- 
tional potential acting on the disk as 



Gnu 



(54) 



where r p denotes the radius vector of the planet location. We 
follow Klahr & Kley (2006) in using a cubic planetary potential 



nu G 



-2 



m p G 



+ 2- 



for d < r sl 



for d > r„. 



(55) 



where r sm = 0.6//. We compute our model using both 
the standard integration method and the proposed FARGO- 
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Table 2. Growth rates of the linear MRI mode. 
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Fig. 11. Profiles of surface density (top) and temperature (bot- 
tom) as a function of radius for the disk-planet interaction prob- 
lem. Solid red and black lines refer to computations obtained 
using FARGO and the standard integration method, respectively. 



MHD scheme (without magnetic field). The employment of the 
FARGO scheme increased the time-step by a factor of about 3.8 
(top panel in Fig.|9]l. 

Owing to our choice of boundary conditions and the pres- 
ence of the planet, we do not fully conserve the total angular mo- 
mentum although the angular momentum fluctuations in time, 
computed with or without orbital advection, are in good agree- 
ment (bottom panel in Fig. [9}. 

In Fig. [10] we show logarithmic maps of density in the equa- 
torial plane after 100 orbits comparing the computations ob- 
tained with and without orbital advection. The differences are 
only marginal. As a more precise test, we also plot the profiles 
of surface density and temperature as a function of radius, aver- 
aged over the azimuthal direction (Fig. 



111. The global slope of 



the profiles as well as the position and width of the gap created 
by the planet are in excellent agreement between both runs. 

3.5. Linear MRI in global disk 

We test the FARGO-MHD scheme on the linear MRI in global 



disks, following the setup presented by Flock et al. ( 2010[ >. This 
test is very sensitive to the numerical consistency and, at the 
same time, the MRI growth rates are affected by the numerical 
dissipation of the scheme. 

An analytical description of the linear stage of MRI was 
given by Balbus & Hawley ( 1991 ), while global simulations as 
well as the nonlinear evolution of the MRI was presented in 



NJmode 


y MI " (std 2 nd ) 


y MKl (FargQ 2 n*) 




16 


0.70 


0.71 


0.74 


8 


0.57 


0.59 


0.69 


4 


0.00 


0.14 


0.14 



Notes. The first column gives the resolution per wavelength, while 
columns 2 to 4 show, respectively, the MRI growth rates obtained using 
the standard scheme, FARGO-MHD with linear reconstruction (Eq.|29|l 
and FARGO-MHD with parabolic reconstruction (PPM, Eq.[3Ql. 



Hawley & Balbus ( 1991 ). The absolute limit of the growth rate 
for ideal MHD is given for the zero radial wave-vectors q r = 
with the normalized wave vector q z 



q, = L Vl6/15V A /Q. 



(56) 



with the Alfven velocity V A = BJ -\j4np. The critical mode 
q, = 0.97 grows exponentially (*F = ^Poe 7 ') with growth rate 
y = 0.75£X For global disk models with disk thickness H, the 
critical wavelength can be rewritten to 



(57) 



with the angular frequency = R and the Alfven velocity 
V A (see also Eq. 2.3 in |Hawley & Balbus|1991 



We use polar cylindrical coordinates (R, (p, Z) with uniform 
resolution in the region < < tt/3, -Rq/2 < z < Rq/2, and 
Ro < R < 47?o, where Rq is the unit length. The initial density 
p and pressure p are constant across the entire disk patch with 
p = l.o, P = c 2 sP IT, c s = 0.1%), and F = 1.00001. The gas is 
initially set up with a Keplerian speed V? = Ro/R. A uniform 
vertical magnetic-field is placed in the region 2Rq < R < 3Rq. 
We choose the strength of the vertical magnetic field to obtain 
the four fastest-growing modes, fitting in the domain at R — 2, 
B z = Bo/n where Bq = 0.055 and n = 4. We use three different 
resolutions [R,<f>,Z] = [224, 168, 64], [1 12, 84, 32], [66,42, 16] 
with a logarithmic increasing grid to measure the capability to 
resolve the MRI wavelength with 16, 8, and 4 grid cells per ver- 
tical scale-height, respectively. We choose V«(z) = Vo&in4z/H 
for initial radial velocity with the vertical size H of the box. 
Boundary conditions are periodic for all variables in the verti- 
cal and azimuthal directions, while zero gradient is imposed on 
all flow quantities at the radial boundaries. 

We employ the HLLD Riemann solver Miyoshi & Kusano 
(2005), piece- wise linear reconstruction and the second-order 
Runge Kutta time-stepping scheme with and without FARGO- 
MHD. The Courant number is set to 0.33 and the spatial recon- 
struction of the electromotive force at the zone edges is carried 
out using the approach described in Gardine r & Stone| ( |2005| l 

In Table [2] we present the growth of radial magnetic en- 
ergy over radius during the linear MRI phase using the stan- 
dard scheme, as well as FARGO-MHD, with piece-wise linear 
(Equation 29 1 and piece-wise parabolic reconstruction (Equation 
30 1. We determine the growth rate from the time derivative of 
the amplitude maxima for Br in Fourier space at each radius. 
Computations obtained with the FARGO-MHD scheme pro- 
vides, at the same resolution, a higher growth rate than the stan- 
dard scheme. This trend has to be attributed to the reduced nu- 
merical dissipation as confirmed by the progressive increase in 
the growth rate when moving from a second-order to a third- 
order interpolation scheme. The observed speedup is around 4. 
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We also show that we need at least 8 grid cells per wavelength to 
correctly resolve the growth rate. The third-order reconstruction 
method in combination with FARGO-MHD reaches with 16 grid 
cells per H, growth rates up to the analytical limit of 0.75. 

3.6. Turbulent accretion disk 

As a final application, we consider a turbulent accretion disk 
in 3D spherical coordinates (r, 0, <p) using the setup BO de- 
scribed by Flock et al. (2011 ). The initial condition consists of 
a vertically stratified structure in a Newtonian central potential 
(5 = — 1 jr with density 



, , / sin (0) - 1 



(58) 



where po = l,R = rsin(0) is the cylindrical radius, and H/R = 
Co = 0.07 defines the ratio of scale-height to radius. The pressure 
follows locally the isothermal equation of state p = c 2 s p where 
the sound speed c s = cq/ y/R. The azimuthal velocity is set to 



Ks= A - 1- 



2.5 

sin(#/ 



(59) 



The initial velocities V r and Vg are set to a white noise perturba- 
tion with amplitude 10~ 4 c s . The simulation uses a pure toroidal 
seed magnetic field with constant plasma B = 2P/B 2 - 25. 

For the sake of comparison, we employ the same compu- 
tational domain, numerical resolution, and boundary conditions 
adopted by [Flock et al. ( 201 l| l. We thus have 1 < r < 10 (in 
AU), tt/2 - 0.3 < < n/2 + 0.3 (approximately ±4.3 disk scale- 
heights), and < (p < 2n. The grid resolution is uniform with 
N r = 384, N e = 192, and N$ = 768 zones in the three directions. 
Buffer zones extend from 1 AU to 2 AU as well as from 9 AU 
to 10 AU. In the buffer zones, we use a linearly increasing re- 
sistivity (up to 77 = 10~ 3 ) reaching the boundary. This dampens 
magnetic field fluctuations and suppresses the interactions with 
the boundary. Our outflow boundary condition projects the radial 
gradients in density, pressure, and azimuthal velocity into the 
radial boundary, and the vertical gradients in density and pres- 
sure at the boundary. For the first run, we employ the HLLD 
Riemann solver, piece-wise linear reconstruction and 2 nd order 
Runge Kutta time integration. The second run is repeated with 
the FARGO-MHD scheme described in this paper. 

As the evolution proceeds, the initial configuration becomes 
vulnerable to the MRI instability, which quickly leads to a turbu- 
lent behavior. Fig. 12 shows a slice-cut of the magnetic pressure 
and density after 600 orbits. 

A direct comparison between the results obtained with the 
FARGO-MHD scheme and those of IFlock et al. 



vided in the three panels of Fig 13 by plotting re 



(2011 



is pro- 
evant volume - 

and time-averaged quantities. For our analysis, we use the range 
between 3 AU and 8 AU (in r), which is unaffected by the buffer 
zones. Temporal averages are taken between 600 and 1200 inner 
orbits. A detailed description of the measurement can be found 
in Section 2.1 in |Flock et aT1 ( |20TT) . 

The volume-averaged values {ass) of the Shakura and 
Sunyaev parameter defined by 



ass 



jp(f -a) dv 

JpdV 



(60) 



stresses even within the first 200 inner orbits, at the linear phase 
of MRI. While the standard computation shows a time-averaged 
a ss value of 5 ■ 10 3 , the results obtained with FARGO-MHD 
reveal a factor of « 2 increase (ass - 9.6 • 10~ 3 ). 

The spectra of B^im) over azimuthal wavenumber m, plotted 
in the middle panel of Fig. 



13 reveals smaller resolved scales 



are plotted, as a function of time, in the top panel of Fig 13 
The FARGO-MHD run shows a much faster increase in the 



(within a factor of two) of the turbulent magnetic field. 

The unstratified global simulations of Sorathia et al. ( 201 1\ 
suggest that the magnetic tilt angle could be a reliable indica- 
tor of numerically resolved MRI. We measure the tilt angle for 
the magnetic field, which is defined by sin20 B = \B r B,p\/B 2 , at 
4.5 AU and plot the time-averaged vertical profile of 0b in the 
bottom panel of Fig. [13] The tilt angle present the highest val- 
ues in the coronal region. Employment of the proposed orbital- 
advection scheme results in a overall larger tilt angle, of * 10° 
in the midplane, compared to the standard computation, where 
* 8° in the midplane. These results suggest that, at the same res- 
olution, the MRI is more accurately resolved when using orbital 
advection. 

With the given grid resolution, chosen to match the results 
of Fl ock et al.| ( |201 l| l, we observe a time-step speedup of ~ 3.75 
when employing FARGO-MHD. Nonetheless, we note that, ac- 
cording to the guidelines given in 5 2.2 this choice may not be 
an optimal one since the cell aspect ratio is 1 : 0.67 : 1.74 
(Ar : rA0 : rA<p), and larger gains may be obtained by suitably 
re-adjusting the number of points in order to have cells with as- 
pect ratios closer to one. 



4. Summary 

We have presented an orbital advection scheme suitable for 
the numerical simulations of magnetized, differentially rotating 
flows in different systems of coordinates. The algorithm has been 
implemented in the release 4.0 of the PLUTO code for astro- 
physical fluid dynamics (Mi gnone et al.|20O7 2012 1 and shares 
the same ideas as proposed by Masset ( 2000J, which consist of 
decomposing, at each time-step, the total velocity into an az- 
imuthally averaged, mean contribution and a residual term. By 
taking advantage of operator splitting, the two contributions are 
carried out as a regular update of the standard MHD equations 
written in the residual velocity and a linear transport step cor- 
responding to a non-integer shift of zones. During the former, 
any dimensionally unsplit scheme may be employed provided 
that additional source terms are taken into account and correctly 
discretized to preserve conservation of both total angular mo- 
mentum and total energy to machine precision. In both steps, the 
magnetic field is evolved using the constrained transport formal- 
ism to maintain the divergence-free condition. 

This approach yields substantially larger time steps when- 
ever the orbital speed exceeds any other characteristic wave sig- 
nal since the Courant condition depends on the residual velocity 
rather than the total velocity. Numerical tests confirm that the 
proposed orbital advection scheme yields results that are equally 
accurate to and less dissipative than the standard numerical ap- 
proach at a reduced numerical cost. The overall gain is problem- 
dependent and can be optimized by a suitable choice of grid res- 
olution in the specified geometry: for the selected problems, we 
observed speedup factors between 4 and 12. 

Parallel-domain decomposition can be performed in all three 
coordinate directions thus allowing the algorithm to be effi- 
ciently employed on large numbers of processors on modern 
parallel computers. The proposed FARGO-MHD scheme is par- 
ticularly suited to large-scale global disk simulations that have 
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only recently become amenable to numerical computations on 
petascale systems. 
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Appendix A: Equations in Different Systems of 
Coordinates 



In this section, we give explicit expressions for the MHD equa- 
tions (|8]l through (TTJ written in terms of the residual velocity 
v' = v — w for different systems of coordinates. Also, to make 
notations more compact, we use m — pv and m' = pv' to de- 
note the total and residual momentum and £' = -v' x B will be 
adopted as a short-hand notation for the (residual) electric field. 



A. 1. Cartesian coordinates 



In Cartesian coordinates (x,y,z), we assume, without any loss 
of generality, that the bulk orbital motion takes place along the 
y direction, i.e., w = wy. The magnitude of w can be defined 
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The source terms on the right-hand side of the momentum and 
energy equations may be written using the reduced forms given 
in Eq. ( 13 1 and ( 14 1. In this case, one has 



S m'. — 



<9<D 

^ dx 

dy 



(A.2) 



S E ' 



<9<D 

~dz~ 



= -pv' ■ VO - pv' (y • Vw) + B y (B ■ Vw) . 



However, following the guidelines given in 4 2.3. 1 



venient discretization is given by using Eq. (21 1 and (22 1. This 



a more con- 



affects only the source terms appearing in the 3'-component of 
the momentum equation and in the energy equation which can 
now be cast as 



S m' y 

S E > 



<9<D 

-p V • (wpv ) + wV ■ (pv ) 

dy 

-pv' ■ VO + wV • (mV - B^B) + wV • {pwv')+ (A3) 



-V • (pv') - V 



pw 



v' - w yn' y v' - B^B} 




Standard 
FARGO- MHD 



A.1.1. Shearing-box equations 

As a particular case, we briefly review the equations of the 



shearing-box model introduced in 5 3.3 By adopting a non- 



inertial frame that co-rotates with the disk at orbital frequency 
Qq, the momentum and energy equations (|3]and|4|i become, re- 
spectively, 



d(pv) 

dt 
dE 



+ V ■ (pw - BB) + Vp, 



pg s - 2Q. z x pv 



(A.4) 



dt 



+ V-[(E + p,)v-(vB)B] = pvg s , 



-2 



() 

SH 



where g s = Q.^(2qxx - zz) is the tidal expansion of the effective 
gravity while q is the shear parameter (Eq 47 1 and the second 
term in Eq. (|A.4[) represents the Coriolis force. The continuity 



Fig. 13. Top: Time history of the accretion stress in the turbu- 
lent accretion disk problem. Middle: Magnetic field spectra over 
azimuth. The FARGO method resolves much smaller turbulent and induction equations retain the same form as the original sys 
scales using the same resolution. Bottom:Magnetic tilt angle tem in E 1- S and VI 
over height. 

either analytically or by averaging in the y direction. In any case, 
w = w(x, z) so that d v w - 0. Writing the equations explicitly 



The derivation of Eq. ( |A.4[ > with FARGO-MHD is done sim- 
ilarly to the previous section with w = -q£loxy leading to 



dm x _ , , „ dp, dm x 
(m x v' - B X B) + ^ + w—± 
at ox oy 

dm' , \ dp, dm' 

^+V-(m'v'-B y B) + ^+w- 



dt 

dm 



dy 
dp, 



dy 
dm- 



+ V -(my -B Z B)+ ^- +w 
dt dz dy 



dE' 



+ V- 



(E'+p t )v'-B(v'-B) 



+ w 



dE' 

dy 



dB x dS z d&y dB x 

dt dy dz dy 

dB y d&' x d&[ d(wB x ) d(wB-) 

dt dz dx dx dz 





S m x 
S m'y 

S HI . 

S E ' 






d(pv') 

dt 
dE' 



+ V ■ (pvV - BB) + Vp, 



(A.5) 



dt 



+ V-[(E' +p t )v' -(v' B)B] = S E >, 



where the source term S,„> and S e< can be shown to be, respec- 
tively, equal to 



(A.l) 



S m , = 2pQ v' v * + pQ (<7 - 2)v x y - pQ^zz 
S E ' = -pv-Q„z - [B y B x - pv'yVxj qtlo ■ 



(A.6) 
(A.7) 



We note that only the vertical component of gravity is included 
in the orbital frame ( Stone & Gardiner 2 010[ ). 
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A.2. Polar coordinates 



A.3. Spherical coordinates 



In polar cylindrical coordinates (R, 0, z), we consider the bulk or- 
bital velocity to be aligned with the azimuthal direction, that is, 
w = CIR()> where Q. = Q.(R, z) is the angular rotation velocity de- 
fined by averaging v$/R along the azimuthal direction. Writing 
the equations in components yields 



dt 1 * * I r d<f> 30 

dm, _ , , „. dp, „dm~ 
— - + V ■ (mj - B,B) + — + Q — - 

dt y ~ - ' dz dcp 





s, 



BE' 
~dt 



(E' +p,y -B(v' B) 



,3/r 

30 



dB R 1 d&' 7 d&j, dB R 

+ — 7T- + — 

8t R d(p dz 30 

dB# dS' R dS' z d(wB R ) d(wB z ) 
~~dt~ + ~dz~ ~ ~dR dR dz~ 

8B Z 1 diRGf) 1 dS' R 8B Z 

+ tt- 1 ^ + ^^r~ 

dt R dR R d<p dif> 



S E > 


0. 



(A.8) 



where, besides the usual divergence operator V • ( ), we have in- 
troduced the "augmented" divergence 



y r F= I d(R F R ) 1 3f> 3F Z 
fi 2 37? 50 3z 



(A.9) 



which avoids the curvature source terms in the 0-component of 
the residual momentum equation, and is more appropriate for 
ensuring total angular-momentum conservation. 

The source terms in the momentum and energy equations can 
be written, using the reduced form Eq. ( p"3) and (14) , as 

30 pvj ' 



- p dR + 



R 



_ p30 , pv R w 
* R d(p R 



In spherical coordinates (r, 0, 0) the direction of orbital motion 
coincides with the azimuthal direction, i.e., w = Qr sin where 
Q = Q(r, 0) is total angular rotation velocity computed by av- 
eraging vJirsmff) in the direction. Writing equations ([8]l 
through ([TT) in components yields 



dm r „ . , ^ dp, dm r 

of or 00 

3m e , 13/?, Jm e 

— + v.(„ w - B „ B)+ -_ +!i _ 



~37 

3B r 



rsinfl 30 



= 

- ? 

= 5 m' 



+ V 



(E' + p,)v' -B(v' -B) 



a 



1 o(sinfl£ ) 
37 rsin# d~9 



1 d6' 
rsinO d<p 
dB e 



30 
dB r 



(A. 12) 



+ o 



S E - 




dt rsinfl 30 r dr d(f> 
dB 1 g(rg^) _ 1 dS' r _ 1 d(rwB r ) _ 1 3(w5 e ) _ 

dt r dr r d6 r dr r d6 
where, in analogy with the polar system of coordinates, we have 
introduced the "augmented" divergence operator 



1 d 



1 



8 



r 3 3r r sin 2 6»d# 

1 



(sin 2 0F 9 ) 



(A. 13) 



rsin# 30 

which is more convenient for expressing conservation of total 
angular momentum. 

The source terms in the momentum and energy equations can 
be written, using the reduced form Eq. ( 13 i and ( 14 1, as 

„2_ B 2 n „2 _ B 2 



dO pv l e - Bg pv\ 

-p— + + 

or r r 



S m z — P 



do 

~dz 



(A. 10) 



S E , = -pv' ■ VO - (pv^v - B^B) ■ Vw + 



PVRV^W - BrB^W 

R 



S E , = -pv' ■ V<I) + 



2 -B\ 



p dip pV B V r - BgB r P^ 

— — h cot 8 

r d6 r r 

P 30 pw 

— ~n in ~ P v ' w ( v '" + cot 6v ^ 

r sin d(f> r 

(pVrVj, - B r B$)w 



(A. 14) 



We note that the total velocity appears in the radial momen- 
tum source term. To restore full conservation of total angular 



momentum and total energy, the formulation given by Eq. (21 



and (22 1 is more appropriate. This leads to 



{pvex'4, - BgB^)w , \ 
+ cot# [pv^v - B^Bj ■ Vw . 

Again, in order to enforce conservation of both total angular mo- 
mentum and total energy, the formulation given by Eq. (21 1 and 



S K = V R ■ (wpv') + wV ■ (pv') 

* R 00 

S E , = -pv' ■ VO + wV R ■ (m'v' - B^B) + wV R ■ (pwv')+ (A.l 1) 



( |22) is the more appropriate to adopt. In this case, one rewrites 
:nergy source terms as 

V • (wpv') + wV • (pv') 



the 0-momentum and energy source terms as 
P 30 



* r sin 6 30 
S E , = -pv' ■ VO + wV • (m'pv' - B,pB) + wV ■ (pwv')+ (A. 15) 



w 



-— V-(pv')-V- 



pw 



v' - w (m'fV' - 



-V • (pv') - V • 



pw 



v' -w (m'^v' - B^fi) 
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